version 13
set more off
global maindir "/home/`c(username)'/Dropbox/Research/Nazis"

use "$maindir/Data/master.dta", clear
eststo clear

* required deps
// ssc install xtsemipar
// ssc install bspline

* older original version of xtsemipar and bspline
* (newer works as well but slight deviation in results)
// net install st0296.pkg
// net install sg151.pkg
// net install sg151_2.pkg

******************** SEMIPARAMETRIC ANALYSIS **********************************

* use only post 1929 data (missings for some outcomes otherwise)
keep if jahr > 1929
keep if id != .

* keep only relevant variables to satisfy maxvar limit
keep infoaratel1lebgeb subdocmel1psh emigration jahr id

* get control function
cap drop vc vcsq vccu
xi: xtreg subdocmel1psh emigration i.jahr, fe vce(cluster id)
* vc linear, control function (squared/cubic in robustness checks)
predict vc, e

* estimate the partialled out residuals, then do nonparametric regression on remainder
capture noisily drop nfit wolfit
xi: xtsemipar infoaratel1lebgeb vc i.jahr, nonpar(subdocmel1psh) generate(nfit wolfit) cluster(id) nograph

* repeat lpoly to get the same parameters
lpoly wolfit subdocmel1psh, ci degree(4)
local b = r(bwidth)
local deg = r(degree)

* calculate marginal effects
capture noisily drop ok
capture noisily drop mfx
bysort subdocmel1psh: generate ok = (_n==1)
dydx nfit subdocmel1psh if ok == 1, gen(mfx)
bysort subdocmel1psh: replace mfx = mfx[1]

* 10% bandwidth variation
local bplus  = `b'*1.1
local bminus = `b'*0.9

capture noisily drop mplus mminus
lpoly wolfit subdocmel1psh, bw(`bplus')  degree(`deg') gen(mplus)  at(subdocmel1psh) nograph
lpoly wolfit subdocmel1psh, bw(`bminus') degree(`deg') gen(mminus) at(subdocmel1psh) nograph

generate epsilon = wolfit-mminus
generate ystar   = .

local niter = 10000

* calculate marginal effects
forvalues i=1(1)`niter' {
    capture drop mstar
    capture drop ok
    quietly replace ystar = mplus + epsilon*sign(invnorm(uniform()))
    quietly lpoly ystar subdocmel1psh, bw(`b') degree(`deg') gen(mstar) at(subdocmel1psh) nograph
    bysort subdocmel1psh: generate ok = (_n==1)
    quietly dydx mstar subdocmel1psh if ok == 1, gen(mfx`i')
    quietly bysort subdocmel1psh: replace mfx`i' = mfx`i'[1]
}

* standard deviation, for each value of x
egen rsd = rowsd(mfx1-mfx`niter')

* normal CIs
generate cilow95 = mfx-1.96*rsd
generate ciup95  = mfx+1.96*rsd
generate cilow90 = mfx-1.645*rsd
generate ciup90  = mfx+1.645*rsd

* empirical CIs
forvalues i=1(1)`niter' {
    gen t`i'=(mfx`i'-mfx)/rsd
}
egen tlow95 = rowpctile(t1-t`niter'), p(97.5)
egen tup95  = rowpctile(t1-t`niter'), p(2.5)
egen tlow90 = rowpctile(t1-t`niter'), p(95.0)
egen tup90  = rowpctile(t1-t`niter'), p(5.0)

capture drop tcilow* tciup*
generate tcilow90 = mfx - tlow90*rsd
generate tciup90  = mfx - tup90*rsd
generate tcilow95 = mfx - tlow95*rsd
generate tciup95  = mfx - tup95*rsd

* simpler and faster, similar coverage when bootstrapping everything using wild
* cluster bootstrap

save "$maindir/Results/semipar-revised.dta", replace


//=================================== GRAPHS ===================================


global maindir "/home/`c(username)'/Dropbox/Research/Nazis"
use "$maindir/Results/semipar-revised.dta", clear


* generate graphs again, just nicer
* this is the same standard local polynomial regression as above.
* could also use different one and then recalculate marginal effects

* np regression fit w/o scatterplot, clipped
capture noisily drop wolfit2
generate wolfit2 = wolfit
replace  wolfit2 = . if wolfit2 >  40
replace  wolfit2 = . if wolfit2 < -20
capture noisily drop lx ly se ul ll
lpoly wolfit subdocmel1psh, ci noscatter degree(4) gen(lx ly) se(se)
gen ul = ly + 1.96*se
gen ll = ly - 1.96*se
replace ll = -20 if ll < -20
twoway ///
    (rarea ul ll lx, color(gs14)) ///
    (line ly lx, color(edkblue)), ///
    graphregion(color(white)) bgcolor(white) ylab(,nogrid) ///
    title("(a) Nonparametric fit", margin(0 0 3 0)) ///
    xtitle("Physicians (per 1,000 of population)") ///
    ytitle("Conditional mean", height(4)) ///
    legend(off) ///
    note("")
graph save "$maindir/Results/Graphs/npfitscale.gph", replace
graph export "$maindir/Results/Graphs/npfitscale.eps", replace logo(off)
!/usr/bin/epstopdf "$maindir/Results/Graphs/npfitscale.eps" --outfile="$maindir/Results/Graphs/npfitscale.pdf"
writepsfrag "$maindir/Results/Graphs/npfitscale.eps" using "$maindir/Results/Graphs/npfitscale.tex", replace

* marginal effects w/ ci, clipped
capture noisily drop up2 low2
generate up2  = up
replace  up2  =  10 if up  >  10 & up != .
generate low2 = low
replace  low2 = -30 if low < -30 & low != .
twoway (rarea up2 low2 subdocmel1psh, sort lwidth(none) lcolor(gs14) color(gs14)) ///
       (line mfx subdocmel1psh, sort legend(off) lcolor(edkblue)), ///
        graphregion(color(white)) bgcolor(white) ylab(,nogrid) ///
        title("(b) Marginal effect", margin(0 0 3 0)) ///
        xtitle("Physicians (per 1,000 of population)") ///
        ytitle("dy/dx", height(4)) ///
graph save "$maindir/Results/Graphs/mfxciscale.gph", replace
graph export "$maindir/Results/Graphs/mfxciscale.eps", replace logo(off)
!/usr/bin/epstopdf "$maindir/Results/Graphs/mfxciscale.eps" --outfile="$maindir/Results/Graphs/mfxciscale.pdf"
writepsfrag "$maindir/Results/Graphs/mfxciscale.eps" using "$maindir/Results/Graphs/mfxciscale.tex", replace

* histogram
hist subdocmel1psh, width(0.05) fraction /// width(0.1) width(0.01)
    graphregion(color(white)) bgcolor(white) ylab(,nogrid) color(edkblue) ///
    title("(c) Density", margin(0 0 3 0)) ///
    xtitle("Physicians (per 1,000 of population)") ///
    ytitle(, height(4))
graph save "$maindir/Results/Graphs/histsubdoc.gph", replace
graph export "$maindir/Results/Graphs/histsubdoc.eps", replace logo(off)
!/usr/bin/epstopdf "$maindir/Results/Graphs/histsubdoc.eps" --outfile="$maindir/Results/Graphs/histsubdoc.pdf"
writepsfrag "$maindir/Results/Graphs/histsubdoc.eps" using "$maindir/Results/Graphs/histsubdoc.tex", replace

graph combine ///
              "$maindir/Results/Graphs/npfitscale.gph"  ///
              "$maindir/Results/Graphs/mfxciscale.gph" ///
              "$maindir/Results/Graphs/histsubdoc.gph" ///
              , cols(1) graphregion(color(white)) ///
                /*title("Semiparametric analysis")*/ ///
                ysize(10)
graph save "$maindir/Results/Graphs/semparscalemfxciscale.gph", replace
graph export "$maindir/Results/Graphs/semparscalemfxciscale.eps", replace logo(off)
!/usr/bin/epstopdf "$maindir/Results/Graphs/semparscalemfxciscale.eps" --outfile="$maindir/Results/Graphs/semparscalemfxciscale.pdf"
writepsfrag "$maindir/Results/Graphs/semparscalemfxciscale.eps" using "$maindir/Results/Graphs/semparscalemfxciscale.tex", replace


graph combine ///
              "$maindir/Results/Graphs/mfxcivc.gph"  ///
              "$maindir/Results/Graphs/mfxcivcsq.gph" ///
              "$maindir/Results/Graphs/mfxcivccu.gph" ///
              , cols(1) graphregion(color(white)) ///
                /*title("Semiparametric analysis")*/ ///
                ysize(10)
graph save "$maindir/Results/Graphs/mfxcivcpoly.gph", replace
graph export "$maindir/Results/Graphs/mfxcivcpoly.eps", replace logo(off)
!/usr/bin/epstopdf "$maindir/Results/Graphs/mfxcivcpoly.eps" --outfile="$maindir/Results/Graphs/mfxcivcpoly.pdf"
writepsfrag "$maindir/Results/Graphs/mfxcivcpoly.eps" using "$maindir/Results/Graphs/mfxcivcpoly.tex", replace
